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Abstract. We show that to correctly describe the position of the critical line in 
' the Kauffman random boolean networks one must take into account percolation 

£H i phenomena underlying the process of damage spreading. For this reason, since 

I , the issue of percolation transition is much simpler in random undirected networks, 

than in the directed ones, we study the Kauffman model in undirected networks. 
We derive the mean field formula for the critical line in the giant components of 
these networks, and show that the critical line characterizing the whole network 
results from the fact that the ordered behavior of small clusters shields the chaotic 
behavior of the giant component. We also show a possible attitude towards the 
analytical description of the shielding effect. The theoretical derivations given 
in this paper quite tally with numerical simulations done for classical random 
' O ' graphs. 

o 
o 

PACS numbers: 89.75.Hc, 89.75.-k, 64.60.Cn, 05.45.-a 

' 1. Introduction 

Almost 40 years ago Stuart Kauffman proposed random Boolean networks (RBNs) for 
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modelling gene regulatory networks T|. Since then, beside its original purpose, the 
model and its modifications have been applied to many different phenomena like cell 
differentiation [2J, immune response [3J, evolution [3], opinion formation [5], neural 
networks [5] , and even quantum gravity problems [7] . 
\ The original RBNs were represented by a set of N elements, J^t = 

{ai(t), <72 (i), crjv'(i)}) each element <Ti having two possible states: active (1), or 
inactive (0). The value of at was controlled by k other elements of the network, i.e. 



a i (t + l)=f i (a il (t),a i2 (t),...,a ih (t)), (1) 

where * was a fixed parameter. The functions /< were selected so that they have 
returned values 1 and with probabilities respectively equal to p and 1 — p. The 
parameters k and p have determined the dynamics of the system (Kauffman network) , 
and it has been shown that for a given probability p, there exists the critical number 
of inputs [13] 
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below which all perturbations in the initial state of the system die out {frozen phase) , 
and above which a small perturbation in the initial state of the system may propagate 
across the entire network {chaotic phase). 

In fact, the behavior of Kauffman model in the vicinity of the critical line k c (p) has 
become a major concern of scientists interested in gene regulatory networks. The main 
reason for this was the conjecture that living organisms operate in a region between 
order and complete randomness or chaos (the so-called edge of chaos) where both 
complexity and rate of evolution are maximized [8] [9] [10] . The analogous behavior 
has been noticed in Kauffman networks, which in the interesting region described by 
eq. ([2]) show stability, homeostatis, and the ability to cope with minor modifications 
when mutated. The networks are stable as well as flexible in this region. 

Recently, when data from real networks have become available [ITJ [12], a 
quantitative comparison of the edge of chaos in these datasets and RBN models has 
brought an encouraging and promising message that even such simple model may quite 
well mimic characteristics of real systems. 

Since, however, one has noticed that real genetic networks exhibit a wide range of 
connectivities, the recent modifications of the standard RBN take into consideration a 
distribution of nodes' degrees P(k). It has been shown that if the random topology of 
the directed network is homogeneous (i.e. all elements of the network are statistically 
equivalent), then the network topology can be meaningfully characterized by the 
average in-degree (fc), and the transition between frozen and chaotic phase occurs 
for [14] : 

Several authors [171 EH] have provided a general formula for the edge of chaos in 
directed networks characterized by the joint degree distribution P{k, q) 

M = I (4) 

where k and q correspond to in- and out-degrees of the same node, respectively. The 
formula (j4j shows that the position of the critical line depends on the correlations 
between k and q in such networks. It is also easy to show that the previous results 
([2]) and ([3]) immediately follow from ([4]) if one assumes the lack of correlations 
P{k,q)=P ln {k)P out {q). 

Very recently, it has been shown by finite size scaling methods (FSS) that 
the critical connectivity {k)^ ss significantly deviates from the value established by 
the Eq. ([3]), even for large system sizes [T5]. More precisely, one observes that 
(k)f ss < (fc)c- To support the observation the authors recall other studies [2U] which 
suggest that gene regulatory networks appear to be in the ordered regime and reside 
slightly below the phase transition between order and chaos in opposite to the theory 
which proposes the critical line to be an evolutionary attractor. 

In the present paper, we suggest other explanation of the observed discrepancy. 
We show (both analytically and numerically) that the discrepancies are due to the 
percolation phenomena, which become important in the region of small values of the 
parameter (k). 

To understand the complexity of percolation phenomena in directed graphs let us 
recall the structure of such a graph (23] [26]. In general, a directed graph consists of 
a giant weakly connected component (GWCC) and several finite components (FCs). 
In the GWCC every site is reachable from every other, provided that the links are 
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treated as bidirectional. The GWCC is further divided into a giant strongly connected 
component (GSCC), consisting of all sites reachable from each other following directed 
links. All sites reachable from the GSCC are referred to as the giant OUT component, 
and the sites from which the GSCC is reachable are referred to as the giant IN 
component. The GSCC is the intersection of the IN and OUT components. All 
sites in the GWCC, but not in the IN and OUT components, are referred to as the 
tendrils (TDs) (see Fig. [TJ. 




Figure 1. General structure of a directed network above the percolation 
threshold. 

Size of all components listed above has doubtless impact on propagation of 
perturbations in directed RBNs. Moreover, GSCC and GWCC start to form at 
different values of the parameter (k) (see Fig. [2^). Although it has been shown 
pZ31 how to find the relative sizes of the components (for example GWCC appears 
when (kq) > (q)), the problem of how to implement the results to the theory of 
perturbation spreading in RBNs is still far from being solved. To make the first step 
in this direction, and to show the importance of percolation phenomena on dynamics 
of RBN we concentrate on undirected case of the model. Although the original RBNs 
have been defined as directed ones, the study of undirected networks significantly 
reduces complexity of the problem (see Fig. [5J) . 




<k> < k > 

Figure 2. Schematic plot of sizes of network components as a function of average 
node degree in a) directed ER graphs and b) undirected ER graphs. 

To this end, we organize the paper as follows. In the next section, we present 
numerical methodology and finite-size scaling of perturbation spreading in RBNs. In 
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section 3 we derive general relation describing position of the critical line in undirected 
RBNs with arbitrary distribution of connections P(k), in the analogy of the mean-field 
theory for directed RBNs [13] . Comparing the theory with numerical simulations we 
show significant deviations between the both approaches. Then an improved treatment 
including percolation phenomena is presented in section 4. A summary of our findings 
is given in section 5. 



2. Critical line in undirected random graphs - numerical simulations 

In order to find the position of the critical line in RBN one has to examine the 
sensitivity of its dynamics with regard to initial conditions. In numerical studies such 
a sensitivity can be analyzed quite simply. One has to start with two initial states 
J2o = °2(0), a N (0)} and ^ = {<7i(0),<7 2 (0), ...,(7jv(0)}, which are identical 

except for a small number of elements, and observe how the differences between both 
configurations Y2t ano - St change in time. If a system is robust then the studied 
configurations lead to similar long-time behavior, otherwise the differences develop in 
time. A suitable measure for the distance between the configurations is the overlap 
x(t) defined as 

1 N 

^(*) = l-^^h(i)-^(*)l- (5) 

i=i 

Note, that in the limit N — > oo, the overlap becomes the probability for two arbitrary 
but corresponding elements, o~i(t) and CTi(i), to be equal. Moreover, the stationary 
long-time limit of the overlap x = lim^oo x(t) can be treated as the order parameter 
of the system. If x = 1 then the system is insensitive to initial perturbations (frozen 
phase), while for x < 1, the initial perturbations propagate across the entire network 
(chaotic phase). 

For numerical purposes we define the probability D that the system is sensitive 
to perturbations 

D _ 2^x(t=T)<x(Q) 1 ^ 

R 

where R is the number of generated networks, and T is the number of system updates. 
In our simulations we take RN = 10 6 and T = 200. The Fig. presents a typical 
example of D dependence on our control parameter (k) for different network sizes. 
Then, we apply finite-size scaling method [23] to determine how the probability D 
scales with the system size. Around some critical point, we predict that systems of all 
sizes are indistinguishable except for a change of scale. This suggests 

D((k)) = m, (?) 

where 

In Eq. |7[). / is one of the functions shown in the figure 3a, {k) c is the critical point, 
and N 1 /" provides the change of scale. Fig. [3b shows how the probability D depends 
on the parameter <j> with fitted parameters (k) c — 1.45 ± 0.04 and v = 2.2 ± 0.1. 
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Figure 3. Probability D against a) control parameter (fe) and b) rescaled 
parameter </> for P = 0.5. 



The other problem which should be noted here is the observation that (k) c 
depends on the number of initially perturbed nodes. In the Fig. [4] we plot the 
dependence of normalized critical connectivity 

(fc) c (A) - (k) c 



(k) c 



(9) 



against the number A of initially perturbed nodes in the network of N = 1000 
elements. For further calculations we choose A = 0.032 TV, since then the error 
in (k) c is less than the error arising in finite-size scaling. 




Figure 4. Normalized critical connectivity against the number of perturbed 
nodes in networks of N = 1000 elements. Lines are shown only for better visibility 
of the presented dependence. 

In Fig. O using the method described above, we show the numerically obtained 
values of (k) c against the parameter p. For p — 0.5 critical connectivity is minimal, 
i.e. (k) c — 1.45. Please note that the size of the giant component for this connectivity 
is about one half of the whole network. One can expect that a large number of 
isolated nodes and clusters can significantly affect the perturbation spreading rate in 
this regime. Moreover, it has been demonstrated [25], that the giant component is 
correlated in sparse networks. In the following, we will show that a mean field theory 
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which does not take into account these percolation and correlation issues, although 
correct for large values of (A;), deviates from numerical results for (k) close to 1. 

3.0 



2.5 

We 2.0 
1.5 
1.0 

0.2 0.4 0.6 0.8 

P 

Figure 5. Phase diagram for undirected RBN model. Points show results 
obtained by numerical simulations. The line is a solution of eq. (15) . 



3. Damage spreading in undirected Kauffman RBN - a simple approach 

In this section, we derive a mean field formula for the critical line characterizing 
Kauffman boolean model in undirected and uncorrelated random graphs with arbitrary 
degree distributions P(k). To this end, we partially reproduce and generalize a simple 
annealed calculations that have been for the first time carried out by Derrida and 
Pomeau [T5] . The case of random directed networks has been studied by Aldana [T5] , 
and also by Lee and Rieger [T?] . 

Thus, let Xi{k, t) corresponds to the probability that a given element i of degree k 
possesses the same value in both configurations ^ t and J2t °f the considered boolean 
network, i.e. <7j(i) = <7j(i). It occurs either when all the k inputs of o~i(t) are equal 
to respective inputs of Oi{t), or when the function cf. |T]), ascribed to the node 
i returns the same value for these two configurations. The first case happens with 
probability 

X(qi,q 2 , ...,qk,t-l) = x(qi,t - l)x(q 2 ,t - 1) . . .x(q k ,t- 1), (10) 

where x(qj, t — 1) represents probability that in the previous time step (t — 1) the jth 
nearest neighbor of i having degree qj was in the same state in the two considered 
configurations. It is also easy to see that the second case arises with probability 
p 2 + (1 — p) 2 , when at least one of the k inputs of o~i differs from its counterpart in 5j 
giving rise to the same values of o~i and 5^. Such a situation, in turn, happens with 
probability equal to 1 — X(qi, q 2 , ■ . ■ , q k , t — 1). Taking all the above together we find 
that the probability Xi(k,t) that o~i(t) = Ui(t) is given by 

Xi (k, t + 1) = X{q u . . . , t) + (p 2 + (1 - pf)X{ qi , ...,t) 

= l-2p(l-p)(l-X( qi ,q 2 ,...,q k ,t)), (11) 

where qi,q 2 , • ■ ■ ,q k stand for degrees of nodes found in the nearest neighborhood of 
the node i. 

The equation (fTT| describes dynamics of a single node i of degree k. In order to 
study boolean dynamics of the whole network one has to average the equation, first 
over the nearest neighborhood of i, next over the whole network. The first step simply 
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means averaging over the distribution P(q\, q 2 , ■ ■ ■ , Qk/k), which describes probability 
that nearest neighbors of i have degrees respectively equal to q\ , q 2 , • • • , % 



x(k,t- 



■i) = i-2p(i-p) (i- E X ^ 



,t)P{ qi ,.../k) ,(12) 



whereas the second step corresponds to averaging of the last equation over the node 
degree distribution P(k) characterizing the whole network. Note, that we have omitted 
the subscript i at x(k, t + 1) in Eq. (fl2|) . After averaging, x(k, t + 1) refers to the set 
of nodes having the same degree k. 

At the moment, before we proceed with our calculations let us outline structural 
properties of the studied networks. At the beginning let us remind that the assumed 
lack of higher-order correlations (e.g. three-point or four-point correlations) means 
that a given link does not influence other links of the considered nodes i and j. 

It translates to the fact that the conditional probability P(qi, ?2> • • • , factorizes 

P(qi,q 2 , q k /k) = P( qi /k)P(q 2 /k) . . . P(q k /k), (13) 

where P(qj/k) describes probability that a node of degree qj is the nearest neighbor 



of a node having degree k. Given the formulas ([TO]) . (fT3|) and (fT5|) . the equation ([12]) 
further simplifies as follows 

c(k) = l-2p(l-p) I 1- [J2x(q)P(q/k)\ j. (14 1 



where, since we are interested in the stationary (i.e. 
equation, we have omitted dependence on time t. 



for t — * oo) solutions of this 



M(y) 




Figure 6. The map y = M(y) considered in the text. The solid line corresponds 
to the situation when the only stable solution is (kx) = (k), i.e. x(k) = 1 for all 
values of k. The dashed line shows the case when the second solution (kx) < (k) 
appears. 



P(Qi/k) = TTTPfe), (15) 



Now, assuming the lack of two point correlations, i.e. 

(ky 

which causes that the nearest neighborhood of each node is the same (in statistical 
terms), and then multiplying both sides of Eq. (fT4"]) by k, and finally averaging the 



Critical line in undirected Kauffman RBN - the role of percolation 



8 



resulting equation over the node degree distribution P{k), we get the desired mean- 
field equation which describes stationary states of the Kauffman model defined on 
undirected and uncorrelated random networks with arbitrary degree distributions 

where (fcx) = J^ fc kx{k)P(k). 

At the moment, note that the state (fcrr) = (k), which in fact corresponds to 
the set of conditions x(k) = 1 for all nodes' degrees k, is always a solution of the last 
equation, see Fig. [5] Note also, that this solution may be stable or unstable depending 
on properties of the considered map y = M(y), where y — (kx)/(k) ([TO)) . In fact, one 
can show that the solution loses its stability, when another solution (kx) < (k) of this 
equation appears. For the first time it happens when 

Um ^ - 1, (17) 
y->i- dy 

where the limit y — > 1~ is equivalent to (kx) — > (k)~. Substituting (|16p into l|17p we 
get the condition for the phase transition between ordered and chaotic behavior of the 
Kauffman model defined on undirected and uncorrelated random network 

^ = 1 risi 

(k) 2p(l-p)' 1 J 

where (k) and (fc 2 ) stand for the first and the second moment of the degree distribution 
P(k), respectively. In the following we briefly analyze the formula for the critical line 
Q18[) in classical random graphs. The case of scale-free networks P(k) ~ fc~ 7 , for 
which the second moment (k 2 ) of the degree distribution becomes important, has 
been analyzed in [31 j . 

Thus, since in classical random graphs (k 2 ) = (k) 2 + (k), the formula (JT5J) 
simplifies 

<*>« = W=rt- L (19) 
In the Fig. [5]one can see numerical simulations of the Kauffman boolean model defined 
on these graphs as compared with the expression (|19[) . In our previous paper [31) we 
have suggested that the visible discrepancy between numerical calculations and their 
theoretical prediction for (k) — > 1 (i.e. for p — > 0.5) may result from the fact that 
(fc) = 1 corresponds to the percolation threshold in these networks. A simple heuristic 
argument behind this statement was the following: because the size of the largest 
component near (k) = 1 is significantly smaller than the network size (the network 
is divided into several disconnected components) , any perturbation cannot propagate 
across the entire system, and the frozen phase is easier achieved. It means that the 
closer percolation threshold (A:) = 1 we are, the more crumbled network (separated 
pieces of the whole system) we analyze, and the theoretical prediction given by Eq. (fT9|) 
works worse and worse. In fact, comparing the general formula (k 2 ) / (k) = 2 [32] for 
the percolation threshold in arbitrary undirected and uncorrelated random network 
with the general expression for the critical line (|18p . one can show that the arguments 
exposed in relation to classical random graphs should also apply for the whole class 
of the considered networks. 

In the next section we show how to adjust the approach presented in this section 
in order to correctly describe properties of the analyzed systems in the whole range of 
parameters, also in the vicinity of the percolation transition. 
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4. The effect of percolation phenomena on damage spreading 

In the following, in order to correctly address the problem of damage spreading in the 
vicinity of percolation transition, that has been outlined at the end of the previous 
section, we use a few important results on percolation phenomena in the considered 
class of networks. To begin with, we recall these results. As we are going to directly 
(i.e. in the course of numerical simulations) check our derivations in classical random 
graphs, together with general formulas describing behavior of arbitrary undirected 
and uncorrelated random networks we also provide the respective formulas for these 
graphs. 

Thus, as we have already mentioned, random graph with a given node degree 
distribution P(k) does not need to be connected. However, if 

W > % m 

that in classical random graphs translates into 

(k) > 1, (21) 

the giant component GC emerges which gathers a finite fraction of all nodes and links. 
The size of the giant component S, i.e. the probability that an arbitrary node belongs 
to GC, is given by the below formula 

S=l-G (tO, (22) 
where u is the solution of the self-consistency equation 

u = Gi(«), (23) 
and 1 — u 2 is the probability that a link belongs to the giant component. The functions 
Go (u) and G% (u) correspond to generating functions of the node degree distribution 
P(k), and the conditional distribution P(qj/k) (p~5|) . respectively. Since in classical 
random graphs Gq(x) = G\{x) = e^ k '^ x ~ l \ the formula (f22|) for these networks 
significantly simplifies 

S=l-e~ {k > s , (24) 

and the expression for u becomes 

u = l-S. (25) 

The general results on percolation transition in random undirected and 
uncorrelated networks outlined in the previous paragraph are already well-known. 
They have been derived by several authors using different theoretical approaches, see 
e.g. [32, 30] . Recently, however, a new interesting results completing our knowledge 
in this subject have been obtained by Bialas and Oles [25]. The authors have shown 
that the neighboring nodes in the giant connected components are disassortatively 
correlated. They have also derived analytic formulas for the node degree distribution 

P*(fc) = P(fc)ir^, (26) 

and the joint nearest-neighbor degree distribution 

-(M)-(M)(i^)=M(iz^!), (27) 

characterizing the giant component. Let us note, that in the limit u — > 0, when 
the giant component covers the whole network S — > 1, the both distributions 
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P*(k) and P*(k,q) respectively converge to distributions P{k) and P(k,q), which 
characterize random uncorrelated networks. The formulas (f2T)]) and (|2~T|) are crucial 
for the further developments of this paper, as they show that although in average 
the considered networks are uncorrelated, in the vicinity of percolation transition 
their giant components are disassortative (note that we still do not know anything 
about higher-order correlations in GCs). Now, since we know that this type of 
correlations makes different spreading- like phenomena more difficult [21] . we expect 
that disassortativity of the giant component is partially responsible for the discrepancy 
observed in Fig. with the crumbling of the system as a whole being the second 
reason. Below, we show that taking these effects into consideration significantly 
improves theoretical prediction for the critical line in the Kauffman model defined 
on random uncorrelated networks. 

Thus, let us study damage spreading within the giant component of the considered 
networks. Knowing properties of this cluster, we can start our analysis from Eq. (fTi|) . 
which is valid for the general class of networks with two-point correlations. The 
conditional probability P*(q/k) for the giant component can be calculated from the 
standard expression [33] 

where 

<fc)*=]TfcP*(fc) = <fc)i-^, (29) 
k b 

is the average degree characterizing this component. Inserting (f2l)]) and (f2"T|) into 
we get 

'1 - U fe +9- 2 ' 
1 - U k 

where P(q/k) is given by (fl"5|) . The last formula ([30]) can be also written in the 
equivalent form 

(31) 



P*(q/k)=P(q/k)( 1 ), (30) 



{k) (l-uQ)J \ 1 

which turns out to be useful in our further developments. 

Now, let us apply the equation (fl4|) to the giant component 

x*(k) = 1 - 2p(l-p) - (^2x*(q)P*(q/k)^ j . (32) 

Due to the complicated form of the conditional distribution P*(q/k) ([30]) . it is 
impossible to deduce on possible solutions of the equation (I32|) in the same way as we 
have done it for the case of uncorrelated networks. However, substituting ([31]) into 
(1321) we obtain 



\k) = 1 - 2p(l -p)L- K(q)w(q, k)P* (q)^j j 



(33) 



where 

( ri\ — T*ln( ... .... 

1 - ui (k) ' 



K (q)= x *(q)-?—«, (34) 
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w{q,k) = — -j— . (35) 



and 

1 - U i +k - 2 
1 — u 

Next, applying a mean field approximation to Eq. (|33[) 

(K(q)w(q, k))* = ]T K(q)w(q, k)P*(q) (36) 



J2<l)P*(q?j ^>(<Z,fc)P*(g)j =«•«;*(*), 



we get the simplified equation 

x*(k) = 1 - 2p(l -p) (l - ( K *w*(fc)) fc ) , (37) 

which after some algebra, consisting in multiplying both sides of this equation by 
(k/(k))(S/(l — u k )) and then averaging it over P*(fc), further simplifies and becomes 
equivalent to Eq. (fl"6)) 

'£m\. (38) 



k*=M*(k*) = l-2p(l-p) ( 1-^(kW(A;)) 



The equivalence of the two equations (JT6J) and (|38[) is visible when u — > (i.e. 
5 — > 1). Then, the parameter k*, see Eqs. (|3"4")l and (|36p . simplifies as follows 

where the averages (. . .}* and (...) have their standard meaning (in our calculations 
'*' always refers to the giant component). This equivalence, also makes possible a 
similar analytical treatment of Eq. (f3"5)) , as the one performed in the reference case of 
uncorrelated networks, compare Eqs. (TT4")) - (fT8"|) ). 

Thus, in order to find condition for the transition between ordered and chaotic 
phase of the Kauffman model defined in giant components of random uncorrelated 
networks we have to check when the solution k* = 1 (|39|) . corresponding to x*(k) = 1 
for all nodes' degrees, becomes unstable. In fact, it happens when 

dM*(n*) , . 

lim — = 1. 41 

dn* v ' 

From the equation (|38[) it follows that the condition has a very simple form 

(k 2 w*(k) k ) _ 1 

(k) ~ 2p(l-p)' 
where 

w* (k) = y, w ^> k ) p * (g) = E 1 1 !^' 2 p * (g) (43) 

9 9 

is defined in Eq. ([3^)1 . At the moment, let us note that in the limiting case of u — > 0, 
the parameter w*(k) — > 1, and the formula (|42|) simplifies to the previous condition 
(d. 



(42) 




Figure 7. Phase diagram for undirected RBN model in classical random graphs. 
Dotted line is a solution of basic Eq. H19I I. Filled points represent numerical 
simulations made for the whole network (the same data are shown in Fig. 5). 
Open points and dashed line correspond respectively to numerical simulations 
and analytic prediction of Eq. 14211 for the Kauffman model defined in giant 
components only. Solid line is the solution of final eq. Q47I . Gray area emphasizes 
the set of parameters where the chaotic behavior, although present in the giant 
component, is not yet visible in the whole network. 



It is easy to check, that in the simplest case of classical random graphs the 
parameter w* (k) is given by 

1 -u-u k - x +u fe+lt - 1 

w fc = — h — \h — k\ — ■ 44 

(1 — u)[l — u h ) 

Inserting the last formula into (|42p , and then numerically solving the resulting equation 
for (k) we obtain theoretical prediction for the critical line of the Kauffman model in 
giant components of these graphs. In Fig. [7] one can see that numerical simulations 
quite tally with the theoretical prediction of Eq. (|12")) . Given the figure, we would also 
like to take note of two other interesting effects related to the analyzed problem. First, 
the critical line characterizing the giant component significantly differs from the curve 
described by the formula (fT8|) . It is shifted towards the numerically obtained critical 
line characterizing the whole network. The observation is in some sense promising, 
as it partially confirms the main proposition of this paper, which states that the 
percolation transition is responsible for discrepancies observed in Fig. [5J The second 
effect concerns mutual relationship between the behavior of the giant component and 
the behavior of the whole network. Since one knows that the giant component makes 
up a macroscopic part of the network (it grows linearly with the network size N, 
and becomes infinite in the thermodynamic limit N — > oc) one could expect that 
dynamics of the whole network should reflect behavior of the giant component. Thus, 
the question is, why the numerically obtained critical line characterizing the whole 
network differs from the theoretical prediction for the giant component. In other 
words, why, for the set of parameters marked by the light gray area in Fig. [7j the 
chaotic behavior of the giant component is not visible in the whole network. 

To solve the problem stated at the end of the last paragraph, let us briefly recall 
what the numerical simulations of the Kauffman model consist in. Thus, in numerical 
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studies we check how the initial perturbation of the system x(Q) = 1 — A (JS|), where 
A <C 1, develops over time. In general, when the parameter x(t = T) < x(0) we 
identify the system as the chaotic one. On the other hand, when x(0) < x(t — T) < 1 
we treat it as being in the ordered phase. In reality, however, due to the fact that 
in the vicinity of the percolation transition the considered Kauffman networks are 
strongly heterogenous, they consist of the giant component which is escorted by a 
number of small tree-like clusters and isolated nodes, the systems should by treated 
more carefully. 



t=0 t=T 




Figure 8. Schematic plot of spreading of perturbation in giant component (GC), 
finite clusters (FCs) and isolated points (IPs). In GC damage spreads, in FCs it 
shrinks, while in IPs it does not change. 

To better describe the situation let us choose the system parameters from the 
region that is marked by the light gray color in Fig. [Jj Then, we introduce a quantity 
Q, which measures chaoticity in the system as a mean damage size caused by a single 
node perturbation. If Q > (<)0 then mean damage size grows(shrinks) in time. 
Condition SI = will allow us to derive the relation for the critical line in the whole 
network. 

Let us now divide the network into three parts: giant connected component (GC), 
finite clusters (FCs) and isolated points (IPs). The figure [5] shows schematically how 
the single node perturbation evolves in time in these three parts of the network. In 
studied range of parameters the giant component behaves chaotically, i.e. the mean 
damage size is larger than initial perturbation and Qgc > 0. On the other hand, the 
small density of connections in finite tree-like clusters does not allow perturbation to 
spread out and Qfcs < 0. Because the state of isolated nodes does not change in time, 
then Sl/ps = 0. Now, if one perturb randomly a set of nodes in the whole network, 
fraction S of perturbations will be located in GC, fraction (1 — S)(l — P(k = 0)) will 
be located in FCs, and the rest of them, i.e. (1 — S)P(k — 0) will perturb isolated 
nodes. Now one can write the condition for transition from the frozen to the chaotic 
state in the whole network: 



sn GC + (i - s){i - p(k = o))n FCs = o, 



(45) 
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where VL, P{k) and S depend on (fc). This equation shows that the ordered behavior 
of small clusters can shield the chaotic behavior of the giant component. Only when 
chaoticity in GC is sufficiently developed, this shielding effect becomes neglected. 




Figure 9. Chaoticity f2 (solid line) in the viciity of the critical point (k)*. (k)* 
and (fc) * are the average node degree in GC and in FCs respectively. Dashed line 
presents linear approximation of f2. 

Now, expanding f2 into power series at (fc} = (k) c 

n = il + -^((k)-{k) c ), (46) 

where fl = in critical point (cf. Fig. one gets the final equation for the critical 
line: 

s((ky - (ky c ) = (s- i)(i - p(k = o))«fc>' - <fc):), (47) 

where (k)f = (k)u (cf. eq.(25) in [25 ). The numerical solution of this implicit equation 
is presented in Fig. [7] as the solid line. 



5. Conclusions 



This study was done to investigate the properties of undirected KBN model in the 
vicinity of percolation threshold. We derived a mean field formula for the critical 
line characterizing KBN model in undirected and uncorrelated random graphs with 
arbitrary degree distributions. We have shown that the results of classical mean field 
theory differ from these obtained by numerical simulations. We have shown also that, 
to explain the discrepancies one has to take into account the effect of correlations 
between adjoining nodes in the giant connected component as well as the effect of 
shielding by finite size clusters. As one can see the problem is not easy even for 
undirected networks. As we have shown in figure 1 and 2 a directedness of the 
network introduces further complications into calculations. Nevertheless, we think 
that a similar approach can be derived even for that case. We hope that the presented 
work will encourage others to pursue these topics in the near future. 
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